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1 Introduction 



The measurement of the different observables related to the so-called four-fermion pro- 
cesses constitutes an important part of LEP2 physics. The most important group of 
observables is related to VT-pair production and decay. It is not only the total cross sec- 
tion we are interested in, but also the wealth of different spin and angular asymmetries. 
This leads to rather sophisticated observables depending on multidimensional distribu- 
tions that are also often strongly affected by experimental cut-offs. There are nine distinct 
W decay channels and, owing to small statistics of ly-pair production, one is often bound 
to combine data of rather different experimental signatures to obtain statistically mean- 
ingful results. The situation of the forthcoming ZZ physics or so-called single-W^ physics 
is rather similar. In all cases, theoretical predictions, to be comparable with the data, 
have to be provided in the form of four-fermion Monte Carlo simulations. 

Monte Carlo integration over four-body phase spaceQ is a non-trivial problem. Why is 
it so? First of all, there is over a hundred different four-fermion final states allowed within 
the Standard Model. They differ by the flavours of final-state quarks and/or leptons. Each 
of these states, even at the Born level, is described by many, often over a hundred, Feyn- 
man graphs. For distinct four-fermion final states one finds thus diametrically different 
structures of singularities depending on the set of Feynman graphs describing the particu- 
lar process. The singularities, or more precisely sharp peaks, have to be carefully tackled 
by the generators to assure their efficiency. They can affect the distribution in one or 
more variables of the particular parametrization of the phase space. As an example we 
can think of the t-channel singularities 1/t, s-channel ones such as resonances, etc. These 
singularities occur in various parts of the phase space. Also, if they are present in more 
than one variable they can affect one another, either destructively or by increasing the 
degree of singularity of the second one. Finally one particular four-fermion process can 
have many distinct singular sub-processes, which need to be treated individually. 

Of course, there are also some similarities and symmetries amongst processes as well 
as individual Feynman graphs. If properly used, they can somewhat reduce the size of 
the problem, but unfortunately its complexity remains intact. At the same time, if one 
departs from the Standard Model and introduces any new particles or interaction types, 
the size and complexity of the problem grows. 

At the top of the above list of difficulties, there is also a strong cut-off dependence of 
the cross-section. This, in the case of electrons close to the beam directions, lead either to 
instability of the Monte Carlo integration or makes the generation ineffective since many 
events are generated outside the detector acceptance. 

In the present paper we show a working solution to the above phase-space integration 
problem. It is based on the multichannel approach that is used in the KORALW Monte Carlo 
code Our aim was to enable the generation of all four-fermion e+e~ /1/2/3/4 final 

^ In reality the case is even more complicated, because of additional bremsstrahlung photons. In this 
paper, however, we restrict ourselves to pure Born- level phase space and do not discuss any issues related 
to photonic corrections. Let us note, only, that in the general case the Born-level phase space can form 
well-defined building-blocks of the whole algorithm. 
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states of the Standard Model over the full phase space or its sub-space, within reasonable 
CPU time. 

At the time of the LEP2 Workshop several dedicated four-fermion Monte Carlo 
codes were presented. As far as the integration method is concerned, some general strate- 
gies can be identified amongst these codes: the adaptive algorithms based on VEGAS- 
type routines 0-0] , the BASES-SPRING algorithm p,!],^, the importance-sampling tech- 
nique 0,|lO|,|ll[] and the mu/iic/ianne/importance-sampling algorithms. In the latter group, 
one finds two approaches: of the EXCALIBUR code [|T^,|l3], followed by ERATO and of 
the KORALW |],|], the subject of this paper. 

The major advantage of the multichannel approach is, in our opinion, that it gives al- 
most complete control over the presampling process. The branches correspond to physical 
configurations in a transparent way. At the same time, since one can implement many 
branches, the structure of the generator can be very rich. No complicated dividing of 
the multidimensional phase space is necessary and no potentially uncontrolled automated 
optimization in the core of the algorithm needs to be performed (although some auxiliary 
optimization at the top levels can be useful). 

The multichannel approach to few-body phase-space integration has been used already 
in the past. For instance, the TAUOLA ||15|- |T^] package for r decays works in such a 
manner up to five-body decays. The FERMI SV |T^ code generates the NC-type four- 
fermion processes in that way also. As it happens with the Monte Carlo algorithms, it 
may be difficult to point to the exact source of the multichannel concept. Traces of it can 



be found in a number of other papers as well. Let us mention just a few of these: |lT9|-p3 

It is instructive to note that there are two "layers" in the structure of multichannel 
algorithms. The core of the algorithm is determined by the internal structure of the 
channels that contain the process-specific information on the actual generation algorithm. 
However, on the top of separate channels there is another, external layer, responsible for 
the "branching" process itself. Recently p4[ this second step has been explained in detail. 

The main advantage of the algorithm presented here is, in our opinion, the following. 
Out of various, already known as well as new, ideas, we show how to create in a systematic 
way, with simple building blocks, a universal and extendable presampler for four-body final 
states. Despite its simplicity, this presampler is shown to work effectively in an actual 
Monte Carlo code KORALW. The important advantage of the algorithm proposed here is its 
modularity and simplicity: all the branches are built out of two rather simple elementary 
pieces. 

In this paper we present a complete description of the multichannel four-fermion Monte 
Carlo algorithm of the KORALW code. We start with the detailed description of the "external 
layer" - that is the branching process - which is quite universal. Using notation of 
ref . , we show how it works in the existing algorithm of KORALW. Starting from section 
3, we discuss the "inner layer" - the actual structure of individual channels dedicated to 
generating the four-fermion phase space. We end the paper with a short summary. 
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2 "Multichannel" Master Formula 

As a basis of the generation of the four-fermion phase space we take a multibranch type 
of Monte Carlo algorithm. In this section, let us summarize the details of the "external 
layer" of four-fermion phase space. Our goal is to calculate the following integral 



d^{P,q)\M{P,q)\^ (1) 



d^iP^l) = 7T4^(2-r^M^->:^^l:^ll(^if , - = 4, (2) 

where P denotes the sum of beam momenta, denote the four-momenta of final parti- 
cles, and q denotes collectively all qk four-momenta. Note that eq. (0) includes all the 
normalization factors, including the flux factor 1/2P^. 

Our strategy in solving eq. (|I|) is to simplify |Mp until it reaches the analytically 
calculable level. Then, we generate this simple distribution and reweight it to an exact 
\M\\ 

Let us then replace |Mp with the crude distribution V'ciJ; in which we introduce the 
structure of branches 

d^{P,q)^CR = dHP^^)T.P^^^^(P^- (3) 

The summation goes over branches, with Nbr denoting the total number of branches. In 
this equation we introduced two set of functions: ipcR ■ "^^^ ^^^^ represents the 

actual (crude) distributions to be generated over the phase space. It is in these functions 
that all the physical information is located. Their specific form is the main subject of 
this paper and will be analysed in detail in the next sections. Here we confine ourselves 
to the formal structure of the arrangement of branches. J7* is composed of the Jacobians 
of the change of variables from four-momenta to angular and invariant-mass variables (as 
defined in formulae below), which we anticipate in each channel i and of which we want 
to get rid (and reintroduce later by appropriate reweighting). They are isolated in eq. (|) 
for simplicity reasons. Finally, f)i are some positive coefficients to be specified later. We 
change variables (-P, 5') into angles and masses (cos^*, 0*, s*) defined for each branch in 
different Lorentz frames. By cos^',(^\s* we denote collectively all the angles and masses 
of a given branch i. The crude distribution becomes 

Nbr 

^ rf«l>Xcosr,f , sOp.V^^^(cosr,f , sO, (4) 
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where we introduced the exphcit form of the Jacobian J'^{cos6_\ (p\s^) in new variables 

J\cose\(j)\s') = = Wk, (5) 



1 1 3 

c?$^(cos^\f,sO = -^j^^ds\ds\Y[dcose]d<P). (7) 

The mass-hke variables s\ and s\ present in the integration element denote squares of 
the invariant masses of systems consisting of two or three outgoing fermions. They are 
present in the definition of the A-factors of the J as well, but in this case s*^ can denote 
squares of masses of the outgoing fermions as well. 

We have found that it is convenient to enlarge integration domains of the branches. 
This can be done for each branch in a different way. We do it by introducing 6j functions 

^^^(cosf ,f = ^^«(cosr,f ,/)e,(l]0 (8) 

and extending the integration areas fi* — ^ ^ci?- The Gj enforce the exact phase-space 
limits. Now, by setting Bj — > 1, we arrive at the final form of the crude distribution 

Nbr 

c^*Xcosr,f , sO|f,^^p.^^^(cos^\ 0\ s^. (9) 

i 

We assume that eq. (^) is integrable analytically, with the result 

j pid^\cose\f,s^)tP},j,{cose\f,s!)=p,Afl,ji. (10) 

The random choice of the branch to be used in the generation of a particular point in the 
phase space is performed with the help of probabilities Pf. 

k 

Later, following ref. [24|, we generate points according to distributions ipc^ of eq. (|^). The 
last step to be done is to reintroduce all the simplifications in the form of the appropriate 
weights to complete the prescription on how to calculate the integral (|l|) . For each channel 
i, the transition from 'ipcp^ of eq. (H) to ipc^ of eq. @) can be achieved with the help of a 
simple step-like weight: 

= ^ = Q^. (12) 
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The whole distribution of eq. is then generated from ipcji with the help of the same 
weight w°'{i), the argument i being the number of the chosen branch for a particular event. 
The distribution generated this way 



Nbr 



y p^JtfRyja^) = ^ (13) 

^ ■n-Kf'^ Nbr ^ ' 



is, up to a normalization factor, equal to that of eq. (^). The normalization (integral) of 
eq. (I) can be calculated as 

Nbr Nbr 

UCR = Y.P^J<fhR=Y.P^J^hR{^''{^)^ (14) 

i i 

''Nbr \ /Nbr \ Nbr 

Y.p'^^cr] f E^. w)).) = {^ym)J2p^^cR, (15) 

where we have made use of the fact that Pi defines a ratio of the number of events 
generated in channel i to the total number of events. With (x)- we denote the average of 
X within the i-th channel and with (p. ... ) we denote consecutive average with respect to 
the branches: 

Nbr ^ Nbr 

• • • ))^) = E ^ (^(^' •••))«' JToo\p, • • • = 5Z (^(^' •••)).• (16) 

i i 

In practice, the external summation over channels is realized with the help of the Monte 
Carlo generation of the channel number with the probabilities Pi and, for every event, 
only that individual weight of the particular channel is calculated and used. Furthermore, 

we do not need to calculate ( . . . ))^\ as an explicit sum of averages {w{i, . . . ))j, 

calculated and stored for each channel separately, as in eq. (0). Instead we calculate 
only one average over the whole sample of generated events of all channels: 

N, 1 ^ ^ 



i events in i all events 

(17) 



Here the argument i denotes that for the event generated with channel i we calculate a 
weight w{i, . . .). 

Including the fundamental matrix-element weight 

. |M(P,g)P 
w [P, q) = -~ — — — (18) 
^cr{P, q) 



Let us also note that in general, in our paper, we will omit the symbol limjv_ 
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and repeating the above steps we end up with the following formula for the integral (|1|) 

Nbr 

d^P,q)\M{P,q)\' = {w{k,P,q))yj,Afl,^, (19) 

(20) 



|M(P,g)P 

w{k,P,q) = w\k)w\P,q)= 



^CR{P,q) 

Let us concentrate for a moment on the coefficients pi. If we rewrite them in a more 
specific form as pi = pi/M^jj^, with the condition = 1, we find from eq. (|ll|) that 
Pi = Pi, so that the Pi have a nice interpretation as the actual branching probabilities. 



The master equation (19) then evolves into 

d<l>{P,q)\MiP,q)\^ = {w{k,P,q)) 

wik,P,q) = |M(P,g)p 



^Uhn^P^q) 



(21) 



Qki^k). (22) 



Much as in ref. we will now discuss a possible variation of the above algorithm. 
Turning back to eq. (^, we modify it by removing the Jacobians Ji. 



Nbi 



d^P,q)^CR = d^P,q) Y,P^^CR{P^g) 



(23) 



and instead introduce them in eq. 



Nbi 



^d¥{cos^,f,s^)piJ\cose\f,s')^[;R{cose^,f,s'). 



(24) 



On the course to eq. (||) we simplify jT^-ij — > 1 and, as a result, weights w"'{i) get modified: 



w;"(z) 



CR 



Accordingly, eq. ( ^iD becomes 



rf$(P,g)|M(P,g)p 



{k,P,q)) 



(25) 



(26) 



wik,P,q) = |M(P,g)p 



■gp,^^^(P,g) 



CR 



j\P,q)Qu{^u). (27) 



Note that here the Jacobians JT" are not summed over branches but calculated for the 
branch of a given event only. As the actual form of the Jacobians is rather simple, we 
should not expect a significant difference in the performance of the two algorithms. 
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Finally, let us remark that both of our algorithms are in fact of the "second" type in 
the classification of ref. p^. In the first one (eq. (pT])), we attribute zero weight to some 
events already inside branches in the form of trivial step functions G,. In the second 
algorithm, weights QiJ^^ are continuous functions between zero and one. 

In the language of measure theory the difference of the two algorithms (defined respec- 
tively by formulae (^) and (p6D ) lies in the choice of the basic measure on the phase-space 
manifold. In the second, more natural case, the Lorentz-invariant phase-space element 
takes this role. 



3 The Branches 

In the previous section we gave the complete general description of the "external layer" 
of the multichannel algorithm responsible for the alignment of separate branches. In this 
section we will present the construction of the actual branches in the case of four-fermion 
final states that we have developed for the KORALW code 0. 

Let us begin with remark. Each branch -0^^ of eqs. (|), (^, ( P^ is intended to describe 
a certain type of singularities that one encounters in the Feynman graphs, but there is no 
unique definition of the V'cr functions. They are dummy functions that cancel out in the 
final result. On the one hand, their role is to mimic as close as possible the complicated 
structure of singularities of the true matrix element, but on the other hand they should not 
be too complicated. First of all they must be analytically integrable over the generation 
area in order to be able to normalize the generator. Secondly, because of the complexity 
of the problem, it is, in our opinion, very useful for them to have a modular structure. 
This allows us to write the program and its documentation in a compact and transparent 
way, leaving room for easier modifications in the future. Finally, these simplifications 
obviously cannot go too far as the ipl^j^ functions have to be quite universal and flexible 
to accommodate a variety of matrix-element singularities. 

How do we fulfil in practice all these, partly contradictory, requirements? How do we 
construct the actual ■j/'cr functions? 

The crucial assumption that we make is the partial factorization of ipcR with respect to 
the angular and mass variables (we skip the branch index i for the next two sub-sections): 

fe(cos^,0,s) = f {s)g {cos e,(j)_,s). (28) 

We request that g{cos6_,(j),s), upon integrating over angular variables, becomes indepen- 
dent of masses s 

J d cos 9 dcpg (cos 9, (p, s) = const. (29) 

Finally, we assume that masses s will always be generated first, before the angles 9_, cf). 

With all these simplifications, one may worry that we are too restrictive. For example, 
it may happen that, for certain configurations of singularities, it would be natural to 
reverse the order of generation - angles first and then masses, or even, further, that one 
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should use some combination of these variables as more "physical" . This approach would, 
however, betray to some extent the modularity of the program: it is always up to the 
author of an algorithm to decide where to draw the line. Our priority was set on simplicity 
and modularity of the algorithm. We believe we achieved those to a large extent. The 
solution presented here may look trivial and simplistic, but the point is that the efficiency 
of the generation is sufficient for our purposes, and we do not need to pay the price of 
formulas several pages long, which may be necessary e.g. in the algorithm of refs. 0,118 
Now we can proceed to the details of the mass and angular distributions / and g. 



3.1 Mass Distributions 

As explained above we intend to generate Si, S2 variables first, and then the angles. Gener- 
ically there are two possible ways of aligning the two mass variables: a "chain" -like and 
a "fork" -like, see Fig. 1. 




7714 



Fig. la: "Chain" Fig. lb: "Fork" 



For each of these configurations we need to give the exact phase-space limits 
corresponding to eq. (^. Then we define an extension Q — > flcR of these limits that 
enable the analytic integration of the distributions to simple results, allowing numerically 
fast normalization to unity. 

For the two cases defined above we continue independently: 
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S2 



S2 




Si 



Si 



Fig. 2a: "Chain"-type phase space Fig. 2b: "Fork"-type phase space 

• The "chain". 

The exact phase-space hmits Q{Q'^) of eq. (H) in this case are given by: 

V^max > + ^/s^ > ^/s^ + 1712, > "^3 + "^4- (30) 

By dropping masses, the above domain Q'-^ extends to the simpler one, of Fig 



2a: 



> Si > S2 > 0. 



(31) 



Next, we assume that si and S2 distributions are identical and that the function 
/(si, S2) factorizes 

1 



/(■5l,S2) 



-/l(si)/l(s2) 



(32) 



This is a somewhat restrictive assumption, but as there will be practically no re- 
strictions on the form of a single distribution /i(si), we can always use a sufficiently 
general form of it. It costs some efficiency, but this is compensated by the (almost) 
perfect representation of the phase-space limits already at a crude level and by the 
speed of the code. 

With our assumptions, the integral J-" can be easily calculated: 



f{Si,S2)dSidS2 



dsifi{si) / ds2fi{s2) (33) 
Jo 

(34) 

Fi{x) = J fi{s)ds (35) 
and the distribution can be normalized, assuming one knows the Fi function. 
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The "fork". 

Also in this configuration we assume factorization of the function /(si, S2) into the 
one-dimensional functions /i(si) and /2(s2): 

f{Su S2) = ^/ J M)f2{s2). (36) 



The phase-space limits (0(fi^) of eq 

max ; 



(37) 



) are the following: 

^ > mi + 1712, y/s2 > "^3 + "^4- 

This area can be extended to the polygonal ficR domain, or explicitly to the sur- 
face defined by the difference between the two squares, as shown in Fig. 2b. The 
integration is easy: 



•^(■Smax) 



fl{Si)f2{s2)dSidS2 



dSidS2fl{Si)f2{s2) 

^ Smax 1 4_ 

\Fl{Smax) - i^l(O)] [F2(W) " i"2(0)] 



Fi[ 



FoiSr 



f,{s)ds, J = 1,2. 



(38) 

(39) 

(40) 
(41) 

(42) 



Finally we make a choice of a suitable form of fj functions. We use a multi-branch 
form of it, each branch a being associated with a basic mass singularity type: 

Nmas £ /■ \ Nmas 



E 



fjais) 



Jq fjai^X^dx 



E 



1, 



(43) 



/ia(s) 



(s + m|j) 

(s + m%y^ log" {{Smax + m%) / {s + m%)), 
{{s - Mlf + MlTir\ 



(44) 



+ m. 



The Nmas denotes the total number of "mass branches" . The parameter m|, takes a role 
of a regulator, preventing distributions from blowing up to infinity at the phase-space 
limit. It should be stressed that m|j does not play the role of a cut-off, and that the 
complete phase space is covered by the generation. The Mk and Ffc are parameters of the 
fc-th resonance. All the functions presented above are easily integrable into the elementary 
functions. 
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3.2 Angular Distributions 

In the next step we generate the angular variables: cos^j and (pj (the index j = 1,2,3 
denotes the number of angle within the same branch). 

The role of the angular variables is to describe the class of t- and u-channel-type 
singularities of the Feynman graphs. There is a number of different transfers that can be 
built, based on two beams and four external final-state four- vectors. In principle each of 
them can develop a singularity in some of the Feynman graphs for some configurations of 
final-state fermions. 

The general angular function g{cos6_, 0, s) of eq. (pHj) is further factorized with respect 
to the angles of each step of the generation. Unlike the mass distributions case, the 
factorization is not complete - each subsequent angle uses previously generated angles 
and masses (in the form of four- momenta). We write it symbolically as: 

5f (cos 9, 0, s) = gi (cos 6*1 , 0i ; s)g2 (cos 6^2 , 02 ; 6*1 , 0i , s)g3 (cos 6*3 , 03 ; 6^1 , 0i , 6^2 , 02 , s) . (45) 

Much as with the mass distributions, we use branching structure to build individual 
gj functions (we always take flat distributions in 0j variables and thus their generation 
decouples) : 



gj{cosej,< 



']! ■■■1: 



Nang 



gja{cos9j,s) 



J^^dcosO J^^ d(f)gja{cos6, s 



^ANG 



with Nang denoting the number of "angular branches". The basic angular distribution 
functions gja are: 

' 1, 



gja{cos9j,s) 



log" {{-tj^rnax + m\)/{~tj + m\)), 



(47) 



-Uj + m\) , 



{-Uj + m\) ^ log" {{-Uj.max + m\) / {-Uj + m^)), 

where tj and Uj are functions of cos 6j and s- variables. The 0j variables are independent 
and generated with a flat distribution. The m\ is, as before, a small mass introduced to 
regularize the distributions at singularity points. 



3.3 Event Construction 

Now, we can continue step by step with the generation and construction of the explicit 
form of the event with the help of the variables generated as explained in the previous 
sections. 
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Let us first define the relation between angles and tj and Uj transfers. For any two- 
to-two sub- process of our generation (see Fig. 1), the transfers tj and Uj are defined 
as 

tj = {Pi-qi,r, «, = (Pi-g2,f, (48) 

where Pi, P2 denote four-momena of the first and second beams, and qi ., q2 denote the 
four-momenta of the first and second outgoing objects. These outgoing objects can either 
be the final-state particles or groups of particles ("intermediate states"). Later on, we 
will use the masses of the objects to label them. Thus [ml], [si] will denote respectively 
the first final-state fermion and the system consisting of second, third and fourth particles 
(case of fig. la) or the system consisting of first and second particles (case of fig. lb). In 
every case, the invariant mass of the object is known, either from the input parameters 
or from the mass generation in the previous step of the algorithm. In this way one has 
an access to all possible transfers tj and Uj. 

Angles corresponding to the invariants are easiest to generate in the local CMS frames 
of outgoing four- vectors {qi., q2^). In such a frame the tj{cos9j, s) and Uj{cos9j, s) can be 
expressed in a standard way (we drop the index j here): 

t = {P^-q^f^ql + P^-2q^,P^ + 2\q,\\P,\cos9, (49) 
u = (Pi-gs)' = g? + P|"2g°P2°-2|gi||Pi|cos^, (50) 

with 





(si2 + ql qi) , 


(51) 


p° = 




(52) 




-S12X {l,ql/si2,qi/si2) , 


(53) 


\P1? = 


^Si2A(l,PiVsi2,P|/Sl2), 


(54) 


S12 = 


{Pi + P2f. 


(55) 



Finally we need to define a series of frames (i.e. four-momenta of group of particles) 
used to construct subsequent angles. This can be, as a matter of fact, decoded from 
figs, la and lb. 

For the configuration of fig. la we start by generating the angle /([m^]. Pi) in the rest 
frame of the two beams Pi and P2. Next, we subtract the generated final four- momentum 
[mi] from the second beam P2 and repeat the previous step for the angle Z([m|], Pi) in the 
rest frame of the pair ([m|], [^2]). Finally, in the third step, after subtracting again [m^ 
from P2 — [^i], we build the angle Z ([^3], Pi) in the rest frame of the pair ([m|], [^4]). 

In the case of fig. lb, the first and second angles, Z ([si]. Pi) in the rest frame of the 
pair ([si], [^2]) and Z([mi],Pi) in the rest frame of the pair ([mi], [m2]), are generated 
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as before. In the definition of the third angle, Z ([m3],P2) in the rest frame of the pair 
Z ([m|], [ml]), the direction of the second beam P2 is chosen as the z-axis instead of the 
first one. 



3.4 Normalization 

The last step is to check out the normalization AA^j? '^^ ^ given ipQji used in the definition 
of the branch i of eq. (p!OD. We need to calculate the integral / d^^ipcR- 

Af'cR = j d^\cose\f,s)ij'cji{cos9\f,s) (56) 
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= / dsids2f{si, S2) / JJ^ d cos 6'j(i0j(7j (cos 6'i, 0i; 6'j, = 1. 
^ ^ j=i 

The normalization to unity follows from the specific form of the / and g functions of our 
choice. The property that normalization factors M^f^ = 1, although not necessary for the 
correctness of the algorithm, simplifies the practical use. Let us elaborate more on this 
point. Thanks to the normalization, and if in formulae ( ^41) and (|47|) only single branches 
are present, the product of functions / and g forms the Jacobian of the transformation 
from the domain to the eight- dimensional unit cube. The inverse transformation 
induces the measure on O^^, which is the (crude) probability distribution. In the multi- 
channel case, when coefficients aj^ and bj^ can all be non-zero, the measure (probability 
distribution) on fi^^ is defined as the weighted^ sum of measures defined by individual 
transformations. 



3.5 Optimization 

Finally, we have to discuss the choice of the coefficients Pi, aj^ and bja- It is these pa- 
rameters that assure flexibility of the algorithm. The coefficients Pi, aja and bj^ are free 
parameters of the generator. Their choice depends on the final state chosen as well as 
on the external cuts used in the calculation of cross-section. An optimal choice of these 
parameters is a delicate and important part of the algorithm. We did it "manually", by 
looking at the physical structures they approximate. Our strategy was to eliminate the 
most overweighted events by iterative adjustments of these coefficients and, if needed, 
adding more branches to the master sum of eq. (^). One can easily imagine other, more 
sophisticated, optimization procedures, for example by assuming functional dependence 
of these coefficients on the type of the final fermions or by allowing for functional de- 
pendence on external cut-offs. This latter optimization may be especially interesting. In 
the case of strongly peaked cross-sections, strong cut-offs can downgrade the efficiency of 
the algorithms by causing excessive rejection, unless compensated by retuning of these 

More precisely, weighted with appropriate products of coefRcients/probabihties aj^ and bj^. 



13 



internal parameters. Yet another option would be a numerical optimization of some sort. 
One can mention in this context the approach of Ref. for example. 



To summarize, such a large number of free parameters makes it more complicated to 
find the optimal configuration of the coefficients, but on the other hand it gives a direct 
access to virtually every singular configuration and gives us a possibility for its direct 
modelling. 

We also have to remember that pj, aja and hj^ are dummy parameters. It means that, 
by varying them, one gets a strong tool for testing the correctness of the integration. 

4 The Algorithm 

At this point we are ready to put together all the pieces and write down the complete 
formula for the overall weight and normalization for our algorithm. To be able to calculate 
a of eq. (|l]), the following steps have to be completed: 

1. Generate the branch number i according to probabilities pk- 

2. Generate point (cos^*,0*,s) according to the ipc^ distribution of the chosen z-th 
branch. 

3. Construct the weight w°'{i). 

4. If w'^{i) 7^ construct the phase-space point (g) out of (cos^*, 0*, s). 

5. Calculate the total weight w{i,P,q) of eq. (^If ) (or (^)) in terms of four-vectors. 



6. Calculate the cross-section as an average of the weight w{i, P, q) 

a={w{t,P,q)). (57) 



A few comments are in order here. 

While constructing a set of all branches, one has to include all the necessary permu- 
tations of external momenta of beams and outgoing fermions. The generic two branches 
— "chain" and "fork", need to be applied to all configurations of external momenta to 
describe all the singularities properly. All together, this leads to over fifty branches in the 
case of four-fermion final states^. 

The algorithm as described in this note is the one with variable weights. By means of 
the standard rejection technique it can provide unweighted events as well. 

^ If one wanted to be even more precise, and counted separately each "basic" branch of eqs. (Q) and 
(|47|), the total number of branches would easily exceed 10^! 
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5 Summary 



In this paper we gave a complete description of the Monte Carlo algorithm for the gener- 
ation of the four-fermion phase space, as needed for LEP2 applications. It is based on the 



"multichannel" approach, as explained elsewhere e.g. in |24|. Strong emphasis is given 
to the construction of separate branches with the help of rather simple modular building 
blocks. This approach allowed for transparent writing of the code and testing, as well as 
for easy extensions and modifications in the future. 

The major disadvantage of the present version, or rather the important place for 
future improvements, is optimization of internal coefficients. In particular, in the case 
of big changes of external cut-offs or centre-of-mass energy, the coefficients may need 
retuning. At present this must be done by hand. 

The algorithm has been successfully implemented in the KORALW Monte Carlo pro- 
gram. Generation can cover the full phase space, or any sub-region defined by cut-offs0. 
This is the case of all four-fermion final states at LEP2 centre-of-mass energies. Numeri- 
cally stable results, with a statistical precision of few per mille, can be obtained from the 
runs easily attainable in a few hours of CPU time of any modern PC. 
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